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1. Introduction 

Relativistic hydrodynamics (RHD) is a basic building block in current efforts in 
numerical relativistic astrophysics and relativity Those programmes provide much 
needed theoretical support to observational efforts focusing on extreme astrophysical 
systems. The development of numerical RHD has followed mostly on the steps 
of numerical non-relativistic hydrodynamics, whose technological significance has 
prompted the generation over the past decades of a large amount of mathematical and 
algorithmic known-how. This knowledge has been historically transferred effectively to 
the relativistic case, the prime example being artificial viscosity techniques g], which 
have been adopted and advocated by the pioneering work of Wilson || . 

More recent techniques are based on deeper mathematical understanding of non- 
linear conservation laws 0, |5|, |(|. The relativistic Euler equations are also readily 
analyzed in conservation form (see [|7| for a modern review of related mathematical 
aspects). Numerical formulations of conservative RHD were first presented in || for 
the one-dimensional case. A multi-dimensional extension, using an explicit "3+1" 
decomposition of spacetime, and valid for a general equation of state (EOS), was 



given in 0. In |T0[ a covariant formulation was presented, adapted to polytropic EOS 
and a specific numerical solution procedure (Roe solver). An alternative formulation, 
restricted to the special relativistic case was given in [|TT] . A wide collection of numerical 
applications based on those approaches is reviewed in [|T2]]. 

Adopting the point of view of non-linear conservation laws liberates the analysis 
from the need to adhere to Newtonian fluid dynamical concepts. This seems to be 
a natural approach when using RHD in general relativistic studies, as is exemplified 
by the ease with which one can extend RHD methods e.g., to lightlike foliations of 
spacetimes - a highly non-Newtonian concept. In |13l we introduced a new covariant 



approach, significantly simpler than |10|, but similarly restricted to perfect fluid EOS. 
This formulation was extensively tested numerically, in spherical symmetry, in particular 
in the context of dynamical black hole accretion ||13|| . The main new results in the 
present paper extend this formulation to a general EOS and introduce a simple linear 
analysis of the RHD equations which allows the tailoring of geometric source terms 
to situations of interest. The discussion unifies and clarifies apparently disconnected 
choices the characterize the literature mentioned above and will hopefully assist the 
future development of the field. 

The paper is organized as follows: In Section ^| we introduce the useful notion 
of linearly equivalent conservation laws and in Section |3| we show how the use of 
the associated transformations allows to tailor the geometrical source structure. In 
particular we develop there a form of the equations suitable for spacetimes with exact 
or approximate Killing symmetries. In Section [| the characteristic structure of the 
standard form of the RHD equations is derived for a general EOS. The inversion of 
the state vector is discussed in some detail. Some classes of EOS which lead to explicit 
inversions are outlined for the algebraically distinct cases of spacelike and null foliations. 
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Numerical applications will be presented elsewhere. Some discussion of related issues is 
included in Section ||. 

Geometrized units (G = c = 1) are used throughout. The metric sign conventions 
follow |TJ|]. Spacetime indices are denoted by small Greek letters and run from zero to 



three. Small Latin indices denote hypersurface coordinates and run from one to three. 
Boldface letters (e.g., F) denote vectors in the fluid state space, which has dimension 
five. Those dimensions are labeled from zero to four, coinciding with the spacetime 
dimension for values up to three. Partial derivatives with respect to a coordinate x are 
denoted as F x and the summation convention is used. 



2. Equivalence classes of conservation laws 

The general form of an N-dimensional system of conservation laws, expressed in a 
coordinate system x M = (x°,x^), where x^ parametrize the hypersurfaces of constant 
time x°, is 

U.. r . • K'.„ S. (1) 

where U is an N- vector describing the state of the system, the flux vectors F J '(x M , U) 
control the time rate of change of the system state within an elementary volume, 
and possible (conservation violating) source terms are denoted by S(a; /1 ,U). The 
independence of the source vector on the state vector derivatives is worth stressing. 

We introduce a iV x N dimensional square matrix G(x M ), which is invertible 
(det(G) 7^ 0). We consider the new state variables tj = GU, satisfying the equations 

%° + F-%3 = S + G ^oG^U + G ^G- 1 ^' , (2) 

where F = GF and S = GS. The linear transformation introduced by G leaves the 
characteristic structure of the original system intact. This can be seen immediately by 
writing the homogeneous version of system ([]]) in quasi-linear form, 

U a o + B'TJ^ =0, (3) 

where 

dF j 

are the Jacobians of the flux vectors with respect to the state vector. The eigenvalues 
of each Jacobian are then determined by the equation 

det{B j - X j l) = , (5) 

where A- 7 denotes an eigenvalue in the direction j and I denotes the unit matrix. The 
corresponding right and left eigenvectors, r 7 and V, are determined by the equations 

(B j - A j I)r J = 0, (6) 
V>(B j =0. (7) 

It follows then from elementary properties of matrices and determinants that the 
transformed system (H) has along the j-th direction the same set of eigenvalues A- 7 as 
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the original system. The new right eigenvectors given by f = Gr, while the new left 
eigenvectors are given by 1 = G _1 l. 

Smooth solutions U(x M ) of the system of equations fll]) obviously lead to smooth 
solutions tJ(x M ) of the transformed system. The same is true for weak solutions, as can 
easily be seen (in the one-dimensional case) from the definition of such solutions || as 
those satisfying the integral relation: 

r+oo r+oo r , r+oo 

/ dx° / dx 1 i!E oU + ^iF 1 = / tfoV(0, ar^U , (g) 

JO J—oo L - 1 J—oo 

where (f){x^) are continuously differentiable test functions of compact support. Pre- 
multiplying this system of integral conservation laws by G immediately shows that weak 
solutions are also transformed properly. Hence from the point of view of a hyperbolic 
set of partial differential equations, all such linearly related systems are equivalent. 

Proceeding to the effects of the linear transformation on the source terms, we 
note that the point-wise properties of G simply reshuffle the source components 
among the equations, but the spacetime dependence of G modifies the RHS of the 
equations in a non-trivial way. Genuine source terms, e.g., such as those introduced 
in Newtonian hydrodynamics to model combustion processes, depend only on the state 
vector and cannot be essentially modified in this manner. In contrast, the regular 
source terms of RHD are highly coordinate dependent and involve a variety of metric 
derivatives. It becomes immediately apparent that an appropriate choice can maximize 
the conservation nature of the system. 

3. Relativistic fluid dynamics 

3.1. Reductions of the equations 

The above discussion is applicable to a variety of relativistic systems but we specialize 
here to RHD. We hence assume a matter current and stress energy tensor corresponding 
to a perfect fluid, i.e., J M = = phu^u 11 +pg fJ,u . In these definitions, is the fluid 

four velocity, which is constrained by the normalization condition u^u^ = — 1, p is the 
rest mass density, p is the pressure, e is the specific internal energy and h — 1 + £ + p/p 
is the specific enthalpy. The pressure is determined by a general two-parameter family 
EOS, p = p(p, e). The relativistic conservation equations in covariant form are given by 

V M T^ = 0, (9) 
V M J^ =0. (10) 

3.1.1. The standard form of RHD Upon introducing a coordinate system, (x°,x l ),one 
obtains a simple form for the continuity equation 

(v^% = 0. (11) 

The standard form of the divergence of the stress-energy tensor can be written as 

(v^n^-r;,^^), (12) 
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where 

r ^A = \g vp {g w ,\ + g\ P ,v - g^x, P ) , (13) 

are the usual Christoffel symbols of the second kind. The set of equations ( p| , p^) 
form a system of non-linear coupled conservation laws which has been the basis for the 
formulations proposed in [JT0|, [IT| , and a "3+1" decomposition of the same equations 
forms the basis of || . We hence call this set of equations the standard form and will 
analyze it further in the sequel. Linear equivalence with many other forms makes our 
analysis applicable to a wider set. 
The state vector is given by 

IF = xf^gT * = (phu°u^ + pg°») , (14) 

u^V^gJ^V^gpu , (is) 

while the flux vectors are given by 

F J > = yf^T^ = 4=g~ (phu^u j + pgM) , (16) 

(17) 

U°) 



F j4 = yf^g J j = yf^g pu J , 
and the geometric source terms are (in this case F j0 



S" 

s 4 



(rft u° + 2rt u fc + rj» f«) 



o. 



(18) 
(19) 



We note that in RHD the source terms depend on the coordinates, the state vector 
and the flux vector. The dependence of the latter on the state vector is implicit, as will 
be elaborated on later. 



3.1.2. The Killing form of the RHD equations It is well known that whereas the 
continuity equation (pi]) is a true conservation law, the divergence of the stress- 
energy tensor will lead in general to true conservation only when a spacetime 
symmetry is present |I5|| . Spacetimes with exact symmetries (e.g., the Kerr spacetime) 
or approximate ones (e.g., quasistationary binary configurations with helicoidal 
symmetry [plj) are of wide interest in fluid-dynamical studies. It is hence of some 
importance to identify ways to maximize conservation in this context. 

A spacetime symmetry is captured by a Killing vector £„ which satisfies the Killing 
equation + = 0. Using all existing symmetries, and complementing any 

missing vectors by appropriate coordinate basis vectors, we introduce four linearly 
independent vectors with a running from zero to three. Upon defining the 

contracted vectors Kfy = & a \T}}, one obtains instead of (Q), 

V,K[ a) = T^V,C ia) , (20) 

which becomes a true conservation law V^lf^s = 0, for each £Yq that satisfies the 
Killing equation. If is an approximate Killing vector in some direction, then the 
corresponding current would be approximately conserved V^K^ « and the related 
source term would capture the, possibly very small, deviations from pure conservation. 
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The state vector U 

U (a) = V 1 

U 4 =V=gJ° 

with the flux vectors given by 



U( a ), U4 in this case is given by 



-g pu° , 



(«) 



-9 C(a) TJ u = \^g {phu 3 i\ a) U v + p&A 



F\ = yf=g~ J j = J=g~ pu j , 

and the geometric source terms are 

S( a ) = ^V^) , 
S 4 =0. 



(21) 
(22) 

(23) 
(24) 

(25) 
(26) 



The linear relation between the standard form of the equations and one adapted to 
Killing symmetries is given simply by the matrix 



G 



£i/(a) 

1 



(27) 



3.1.3. Optimal source terms for RHD in general spacetimes The considerations in the 
last paragraph lead naturally to the case also important for numerical relativity, namely 
the case of no symmetries. Upon introducing a tetrad adapted to the coordinates, and 



using Kfc 



where 



e v {a) T^ the form © results to 



-A p ( 



f —gT P 



pi 1 



1 



g»\g x \v 



With the definition of the state vector according to 



the flux vectors 

F J , = < 
Y\ = , 

and the source terms 

S, 





r —g~T° lJi ■■ 
^g-J° = 

f=gT j » 

r~gp = 



g (phu°u^ + p8l) 




-5 pu 



g (phu J u fl + p5 J ^ , 

j 



-g pu 



So 



(A°,U + A*,U fc 



ALFJ + A* Ft) 



(28) 
(29) 

(30) 
(31) 

(32) 
(33) 

(34) 
(35) 



we obtain a second form for the equations of RHD. Note that for notational economy 
we have used a lowercase state vector index in the above expressions. With the explicit 
substitution of the perfect fluid stress-energy tensor, we obtain 



phu x u p g Xp ^ - pi^f^g ), M 



(36) 
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The source terms in this case involve considerably fewer summations of metric 
derivatives compared to the standard form (0). The mixed form V^T^ = has 
been used before as the starting point for non-conservative approaches to the RHD 
equations ||. We point out here that the the related source simplification should benefit 
conservative formulations as well. The relation of the two forms is of course a lowering 
of the free spacetime index, which is captured by a linear transformation of the form 





G 



9 ni> 




(37) 



We also include here the source terms in the notation of the "3+1" spacetime 
decomposition 
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which is a commonly adopted starting point for developing 
algorithms for dynamical evolutions of spacetimes. Assuming a foliation with spacelike 
surfaces having a unit normal n M , the four-metric is decomposed as 

9nu = liiv - ?fyn„ , (38) 

where 7^ is the 3-metric of the hypersurfaces. Upon introducing the lapse function 
N and the spacelike shift vector N l |15|], the evolution proceeds along the vector field 
ft* = Nn^ 1 + N 11 and the components of the metric read explicitly 

2oo= -N 2 + ll3 WN^ (39) 

g 0j = ltJ N\ (40) 

9ij = lij • (41) 

The state vector we use differs from the usual one (see, e.g., || |17|) m that it 
does not explicitly use the above decomposition. For example, U° = T 00 in contrast to 
p H = T^ u n^n u . The state vector variables hence lose their usual meaning as observables 
in the instantaneous Eulerian rest frame. In a generic spacetime this frame is of no 
special significance and simply reflects the particular choice of lapse and shift vector. In 
compensation, our choice of state vector is valid in the absence of a spacelike foliation. 

The volume element now reduces to \f—g = Ny/^y, where 7 is the determinant of 
the 3-metric, and the source terms for equation ( p8|) have the following explicit form: 

0\2 



^N 2 N^ u y P h 



u 



0\2 



+ N^u°ph ltJ N\ u 



u°N 3 + u 3 



(42) 



4. Further analysis of the RHD equations 

In contrast to non-relativistic hydrodynamics, the relativistic theory exhibits a non- 
linear algebraic coupling of all equations (including the continuity equation) through 
the velocity normalization condition. This feature has wide ranging implications for the 
structure of the theory, and in particular: 
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• The computation of fluxes from the state vector. In non-relativistic hydrodynamics, 
the fluxes can be written as explicit functions of the state vector. In RHD both 
fluxes and state vector are more properly seen as algebraic functions of suitable 
"primitive" variables, as is obvious e.g., from equations ( |I"4[|T7j ). 

• The analysis of the local characteristic structure. The implicit dependence of the 
fluxes on the state vector in RHD requires, in turn, that the analysis of the Jacobians 
of the fluxes with respect to the state vector be done through a set of intermediate 
variables. 

A unique set of such intermediate variables employed in both procedures is 
commonly used in the literature without explicit statement. There are few constraints 
on the choice of primitive variables w. An appropriate choice influences the algebraic 
difficulty of analyzing the characteristic structure of the system, a key ingredient of state- 
of-the-art Riemann solver based numerical schemes for non-linear conservation laws HI. 



The rest mass density p is a common candidate in all proposed sets |]TT] , |9j. With our 
choice of representation (e.g., equations ( |i"2| ) or (^)), the most appropriate velocity 
variable is u l , although a lower index would lead to simpler calculations in the latter 
case. The set must be completed with the choice of an additional thermodynamical 



variable. Both the specific internal energy |J and pressure |IjJ have been used. The 
choice of pressure leads to slightly simpler analysis. A reasonable compromise between 
simplicity and maintaining continuity with past work is hence the choice w = (p,u l ,e). 

4-1- Characteristic structure 

Using the intermediate variables w the conservation law ([!]) is rewritten as a quasi-linear 

system 

A°w )X o + A j w iX j = , (43) 



where 



d\J d¥ j 

— , A 3 = —— . (44) 
dw ' dw 



Hence, upon introducing the intermediate eigenvalues X 3 

det{A 3 - X j A°) = , (45) 

and the intermediate right and left eigenvectors, 

{A j - X j A°)r 3 = (46) 

lj(A 3 - A J A°) = 0, (47) 

elementary algebra establishes that X 3 = X 3 and the right eigenvectors of the matrix B J 
are given by r J = A°r^ . The left eigenvectors of this matrix are given simply by lj = lj. 
The A^ 1 matrices for the equations of RHD are 

Yu°u° + ng 00 2php iU ° Zu°u° + xg 00 

Yu k u° + Kg k0 ph(5 k u° + u k pi) Zu k u° + xg k0 , 
ppi u° 
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Yu°u j + Kg°i ph(niU j + u°5j) Zu°u j + X9 0j 
Yu k u j + ng kj ph(S k u j + u k Sj) Zu k u j + X9 kj 







ft 



(49) 



where Y = p + k, Z = 1 + e + x an d 
du° 



Hi 



Ui dp dp 

— , X = a- > K = ^~ 

M C7p OS 



(50) 



We choose a coordinate direction, which we label £ 1\ The other two coordinate 
directions are then denoted by A = (2, 3). 
The matrix A 1 — A 1 A has eigenvalues 

Aq = v 1 (triple) , 



and 



where 



with v 1 



Ai 



1 



l-cftl + L) 



-Mc 2 + v 1 {1-c 2 s )±c s Vd 



D = c 2 s (M 2 - LN) + (1 - c 2 s )(N - 2Mt> 1 + L^ 1 ) 2 ) , 
^j. The local sound speed is denoted by c s and satisfies 
he 2 



L P 

p 2 



and the following shorthand notation was also used: 



5 



oo 



it 



0^2 



M 



9 



01 



0^2 



N 



9 



ii 



0^2 



(51) 
(52) 

(53) 
(54) 

(55) 



A complete set of linearly independent right-eigenvectors (r 1 = A ? 1 ), is given by 



r o,i 



ro,A 



a 



i5(i + A 8 ),o 



A 



Ut 



+ u A 



oi 1 ll T 

«°,« 1 ,0,0, r -- 
ri a. 



with the definitions 



(56) 
(57) 
(58) 

(59) 

(v* - \ ± )(1 - <*) - <*(M - \ ± L) ■ (60) 
We note that the t± eigenvectors are unique up to normalization, whereas the 
r o,i? r o,A vectors can be any set spanning the degenerate subspace. For a perfect 
fluid EOS these expressions coincide with the ones reported in [HJ. Note that in 
that restricted case a — 1. The spectral decomposition given above applies to a 
chosen direction j. Since j is arbitrary, to obtain similar expressions for the remaining 



a = 1 + e — —p . 

K 
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directions, it suffices to specialize them accordingly, e.g., obtain the eigenvalues from 
expressions (|51| ) and (|52| ) with substitution of the desired direction, and permutation 
of the corresponding eigenvectors. Complete sets of eigenvectors for other versions of 
the equations are obtained by a straightforward multiplication with the corresponding 
G matrix and are not reproduced here. 



4-2. Inverting the state vector 

The general statement of the recovery of the primitives fields is to find values for e.g., 
w = (p,e,u l ), given a set of conserved variables U. This procedure is part of any 
solution algorithm that uses the conservation form of RHD, i.e., it is not related to the 
characteristic decomposition of the system. 

We focus on the standard form of the equations, in which case U = 
y/—g (T 00 ,T 01 , J°). At first sight, this requires the inversion of the system of six non- 
linear algebraic equations 

T 00 =l(2. + l+e)(J )*+pg 00 i (61a) 
P P 

T 0i = (^ + l + e)f ) u i + P g oi , (616) 

V = 9oo(J°) 2 + 2pg oi J°u i + p 2 g ij u i u j , (61c) 

V =P(P,£)- (61c?) 

We point out that even in the most general case the size of the non-linear system 
can be reduced, with the elimination of the velocity from the unknowns. We introduce 
the tensor S^ a = g vp T^ v T pa ' . Inspection of the S 00 component of this tensor shows 
immediately that it is only a function of conserved variables. Together with T 00 and 
the EOS, we have a reduced system which reads 

S™ = (1 + 1 + e)(E - 1 - e)(J ) 2 + p 2 g 00 , (62a) 
P P 

T™ = ±(P + l + E)(J°) 2 +pg 00 , (626) 
P P 

V =p(p,e). (62c) 

The use of a general EOS which may be only available in tabulated form implies that 
the reduced system ( |62Q| - |t)2c) ) is in general to be inverted numerically with an iteration 
procedure. A description of a typical iterative procedure (applied to equations ( |61g{ - |61q 



in the special relativistic limit) is given in ||18|| , where estimated values for (e,v\p) are 



used to start a non-linear iteration of the system. In [|L0| several other procedures are 
discussed for the special case of a polytropic gas. 

It is of interest to point out that further reduction of the system to a binomial 
equation, and hence with an explicit and convenient solution, is possible for special 
classes of analytic EOS. Assuming that the EOS if of the explicit form p = pF{h) allows 
the further manipulation of ( |62q| - |32d) to obtain 



(T m ) 2 S m - h(2F(h) - h)(J°) 2 (T 00 ) 2 - (S 00 - h{F{h) - h)(J ) 2 ) V° = , (63) 
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which is a non-linear equation for the enthalpy. An explicit relation between conserved 
and primitive variables, which rests on the ability of solving equation fl63|), has an 
impact on the efficiency of the numerical code, as it eliminates an iterative process that 
is required, at least once per each spacetime point. 

It is seen immediately that the metric component g 00 stands out as having special 
significance in this algebraic equation. Indeed, as already pointed out in [H| any 
null foliation (characterized by g 00 = 0) leads to explicit solutions for h in terms of 
conserved quantities, in the case of a polytropic equation of state p = (T — l)pe, where 
T is the adiabatic index. The most general case explicitly reducible to a binomial 
is F(h) = a/h + (3 + jh, where (a, /3, 7) are constants characterising the fluid. The 
polytropic gas is the special case with a = 0, 7 = — (3 — (T — 1)/T. 

For spacelike foliations (g 00 7^ 0), the situation is slightly different. The choice 
F(h) = a/h + (3 + h leads to a binomial for h, whereas the choice F(h) = a/h + 7/1 
leads to a binomial for h 2 . None of those cases includes the polytropic gas, but see the 
arguments of [jnj for an example of the latter case, with 7 = —a = 1/4. 

Once the enthalpy is obtained, the other variables follow straightforwardly, e.g., the 
velocity follows from 

(p + p + pe)J° v ; 

5. Summary and concluding remarks 

In view of the increasing importance of conservative RHD in numerical applications 
in relativistic astrophysics and relativity, we explored the corresponding framework 
in some detail. The algebraic complexity of the systems involved suggests that links 
between specific manifestations of the equations be established at the outset. This is 
accomplished with the introduction of simple linear transformations of the equations, 
which leave the local characteristic structure invariant but have non-trivial impact on 
the source terms. It is pointed out that for spacetimes with exact or approximate 
symmetries, but also in the general case, the equations can be modified to capture the 
conservation property in an optimal way. 

The local characteristic structure of the RHD equations in a general spacetime 
foliation (spacelike or lightlike) is analyzed for a general equation of state (EOS). This 
extends previous results restricted to the polytropic case. Special classes of EOS are 
identified for both spacelike and null foliations, which lead to explicit inversion of the 
state vector and computational savings. In a lightlike foliation, the commonly used 
polytropic EOS is included in the explicitly invertible cases. 

Conservative RHD techniques have been applied only recently to general 3D 
spacetime evolutions [XT]- It seems worthwhile to investigate the benefits of the increased 



control over source term structure presented here, in situations of current interest, e.g., 
the study of neutron star binary coalescence. 

The presented framework is uniquely suitable for strong-field simulations using 
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lightlike foliations of the spacetime. Lightlike foliations attached to the exterior of 
spacelike surfaces have been suggested as an effective way for providing global spacetime 
solutions |20| (for a review see plj). The whole approach has been dubbed Cauchy- 
characteristic matching (CCM). Fluid evolution in the CCM context has already been 



investigated P2[ . The RHD framework proposed here is the natural candidate for 
providing economical and state of the art numerical fluid evolution in both foliations, as 
it is form invariant with respect to the slicing. Algorithmically, the only change between 
spacelike/null domains would be the routine used for the recovery of primitives. 

A recent demonstration of long-term stable numerical evolutions of a single black 
hole based on null coordinates |23| opens the possibility for applications to non-vacuum 



single black hole environments solely within a lightlike framework. A full implementation 
of the present formulation in the case of spherical symmetry has been presented |[1 3|| , 



illustrating the ease of studying black hole growth through accretion using an ingoing 
null foliation. Three dimensional implementations for fixed black hole backgrounds in 
null coordinates are also available and results will be discussed elsewhere. 
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